function [P_new,cy1,cy2,cy3,w1,w2,w3] = Value_Iteration_3state(wgrid,P_old,cy1_old,cy2_old,cy3_old,w1_old,w2_old,w3_old,param)

betta = param(1);
delta = param(2);
gamma = param(3);
pi1 = param(4);
pi2 = param(5);
pi3 = param(6);
e1 = param(7);
e2 = param(8);
e3 = param(9);
ey1 = param(10);
ey2 = param(11);
ey3 = param(12);
eo1 = e1-ey1;
eo2 = e2-ey2;
eo3 = e3-ey3;
v1aut = param(13);
v2aut = param(14);
v3aut = param(15);
wmin = param(16);
wmax = param(17);
nw = param(18);

if nargin < 3, cy1_old=0.1*ones(numel(wgrid),1); cy2_old=0.1*ones(numel(wgrid),1); end

%Preallocation
P_new=zeros(nw,1);
cy1=zeros(nw,1);
cy2=zeros(nw,1);
cy3=zeros(nw,1);
w1=zeros(nw,1);
w2=zeros(nw,1);
w3=zeros(nw,1);

w_thr=max(wgrid(((P_old(1)-P_old<10^(-10)))));

% Bounds for c and w
lb = [0;0;0;w_thr;w_thr;w_thr];
ub = [ey1;ey2;ey3;wmax;wmax;wmax];

    function res = value_unconstrained(c)
            
        if (c(1) > 0) && (c(2) > 0)&& (c(3) > 0) && ((e1-c(1)) > 0) && ((e2-c(2)) > 0)&& ((e3-c(3)) > 0)
    
            res=...
            -(pi1.*(utility(c(1),gamma)+(betta./delta).*utility(e1-c(1),gamma)) ...
            + pi2.*(utility(c(2),gamma)+(betta./delta).*utility(e2-c(2),gamma)) ...
            + pi3.*(utility(c(3),gamma)+(betta./delta).*utility(e3-c(3),gamma))+...
            delta.*(pi1*interp1(wgrid,P_old,c(4),'spline') + pi2.*interp1(wgrid,P_old,c(5),'spline')+ pi3.*interp1(wgrid,P_old,c(6),'spline')));
            
        else 
            res=100000;
        end
   
    end

    function [ineqcon,eqcon] = constraints_fun(c,omega)
        
        eqcon=zeros(1,1);
        ineqcon=zeros(3,1);
       
        eqcon(1)=((v3aut-utility(c(3),gamma))./betta) - c(6);
        ineqcon(1)=(omega-(pi1*utility(e1-c(1),gamma)+pi2*utility(e2-c(2),gamma)+pi3*utility(e3-c(3),gamma)));
        ineqcon(2)=((v2aut-utility(c(2),gamma))./betta) - c(5);
        ineqcon(3)=((v1aut-utility(c(1),gamma))./betta) - c(4);
        
    end

%Speed-up: start from the last maximizer
startguess=[cy1_old cy2_old cy3_old w1_old w2_old w3_old];
optionsset=optimoptions('fmincon','Display','off','TolX',10^(-8),'TolFun',10^(-8),'Algorithm','sqp');

%CONSTRAINED VERSION
for k=1:nw
    guess=startguess(k,:)';
    [out,P_new(k)]=fmincon(@(x)value_unconstrained(x),guess,[],[],[],[],lb,ub,@(x)constraints_fun(x,wgrid(k)),optionsset);
    cy1(k)=out(1);
    cy2(k)=out(2);
    cy3(k)=out(3);
    w1(k)=out(4);
    w2(k)=out(5);
    w3(k)=out(6);
    
end

%Minimizers want to minimize - thus here must be negated
P_new = -P_new;

cy1(numel(wgrid));
cy2(numel(wgrid));
cy3(numel(wgrid));
w1(numel(wgrid));
w2(numel(wgrid));
w3(numel(wgrid));
P_new(numel(wgrid))
end